n-body dynamics of stabilized vector solitons. 
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In this work we study the interactions between stabilized Townes solitons. By means of effective 
Lagrangian methods, we have found that the interactions between these solitons are governed by 
central forces, in a first approximation. In our numerical simulations we describe different types of 
orbits, deflections, trapping and soliton splitting. Splitting phenomena are also described by finite- 
dimensional reduced models. All these effects could be used for potential applications of stabilized 
solitons. 



The cubic Nonlinear Schrodinger Equation 
(NLSE) is among the most important physical 
models in the field of nonlinear waves. Besides its 
fundamental value as a first order nonlinear wave 
equation, it is an integrable model in the one di- 
mensional case [l| and represents many different 
physical systems: from laser wavepackets propa- 
gating in nonlinear materials to matter waves in 
Bose-Einstein Condensates (BEC), gravitational 
models for quantum mechanics, plasma physics or 
wave propagation in geological systems, between 
others [2, 3, 4j. In this paper we consider the vec- 
tor version of the NLSE and study the dynamics 
of some particular solutions which are stabilized 
by means of a periodic modulation of the nonlin- 
earity. 



I. INTRODUCTION 

One of the types of NLS equations arising frequently in 
the apphcations is the two-dimensional cubic Nonlinear 
Schrodinger equation, which is of the form 



du 



1 



u(r,0) = uo(r) G H\m.^) 



(la) 
(lb) 



where u(r, i) : M x IR+ — > C is the complex wave ampli- 
tude, A — jdx^ -I- /dy^ and g{t) is a real function 
(the nonlinear coefficient) so that if g < the nonlin- 
earity is attractive whereas for g > the nonlinearity is 
repulsive. 



When g is a real constant Eq. (|la|l is the cubic NLSE, 
which is one of the most important models of mathemat- 
ical physics. 

It is well known that for 5 < if = |uoP, is above 
a threshold value Nc, solutions of Eq. (|la|) can self- focus 
and become singular in a finite time. This phenomenon 
is called wave collapse or blowup of the wave amplitude. 
More precisely, there is never blowup when N < Nc but 
for any e > 0, there exist solutions with iV = iV^ + e for 
which blowup takes place . 

Eq. Hla() admits stationary solutions of the form 
u{r,t) — e*^*<i>^(r), where ^fi{r) verifies 



0. 



(2) 



As it is precisely stated in when g is negative, for each 
positive II there exists only one solution of Eq. which 
is real, positive and radially symmetric and for which 
/ |<i>^pdr has the minimum value between all of the pos- 
sible solutions of Eq. Moreover, the positivity of fi 
ensures that this solution decays exponentially at infin- 
ity. This solution is called the ground state or Townes 
soliton. We will denote it as Rfj.{r) which satisfies 



AR^ - 2iiRf, - 2gRl = 
lim R^{r) = 0, i?' (0) = 0. 



(3) 
(4) 



From the theory of nonlinear Schrodinger equations it 
is known that the Townes soliton has exactly the critical 
norm for blowup Nc, therefore, it separates in some sense 
the region of collapsing and expanding solutions. More- 
over, the Townes soliton is unstable, i.e. small perturba- 
tions of this solution lead to either expansion of the initial 
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data or blowup in finite time. This instability is an essen- 
tial feature of these type of equations which has its origin 
in the two-dimensional nature of the equation and makes 
it essentially different from its one-dimensional version, 
in which stable solitary- wave solutions exist, the so-called 
solitons. 

From the physical point of view the existence of this 
instability implies that no localized solutions of the 2D 
NLSE exist. This is why there has been a great inter- 
est on the case where g is not constant but a continuous 
periodic function of which has arisen recently in differ- 
ent fields of applications of Eq. Hla|l . The intuitive idea 
is that (oscillating) bound states could be obtained by 
combining cicles of positive and negative g values so that 
after an expansion and contraction regime the solution 
could come back to the initial state. In this way some 
sort of pulsating trapped solution, i.e., a breather^ could 
be obtained. 

This idea was first proposed in the field of nonlinear 
Optics 0. In that context, a spatial modulation of the 
Kerr coefficient (the nonlinearity) of the optical mate- 
rial is used to prevent collapse so that the wavepacket 
becomes collapsing and expanding in alternating regions 
and is stabilized in average d, Q . The same idea has 
been used in the field of matter waves in Refs. Cailll. 
In Ref. ^3 some general results are provided for generic 
forms of g{t). Also in Refs. it has been argued 

that the structure which remains stabilized is a Townes 
soliton. The stabilization of more complex structures in 
the framework of Eq. (|la|l constitutes an open problem 
since other solutions such as those of vortex type cannot 
be stabilized [1IT1|. 



II. THE VECTOR NLS AND STABILIZED 
VECTOR SOLITONS 

In this paper we explore the vector version of Eq. (|la|l . 
which is of the form 

du 1 

i-^ = --V^Uj+g{t) (ajiluip + .. . + aj„|u„p)uj, 

(5) 

where j = l,...,n, uj are the complex amplitudes 
A — d'^/dx'^ + d'^/dy'^, a.jk £ K are the nonhnear couphng 
coefficients and g{t) is a periodic function accounting for 
the modulation of the nonlinearity. 

Eqs. are the natural extension of the Manakov 
system ;15] to two transverse dimensions and an arbi- 
trary number of components. In Optics, for spatial soli- 
tons, t plays the role of the propagation coordinate and 
Uj are n mutually incoherent laser wavepackets. One- 
dimensional Manakov-type models have been extensively 
studied in nonlinear Optics, mainly due to the poten- 
tial applications of Manakov solitons in the design of all- 
optical computing devices j^. In BEG these equations 
(with an additional trapping term) describe the dynam- 
ics of multicomponent two-dimensional condensates, Uj 



being the wavefunctions for each of the atomic species 
involved [THIl. 

Some features of this model have been described in 
Ref. In particular it is clear that if Uj{r,0) — 

Rp,{\\r — rjW), where is a Stabilized Townes Soliton 
(STS), and the centers of the distributions are much more 
separated than the width of the Townes solitons, then we 
may have states with STS on each component. Because 
of the Galilean invariance we can also construct solutions 
propagating with uniform velocity along straight trayec- 
tories Tj (t) . Again, if the trayectories are well separated 
it is reasonable to expect that (as it happens in the case 
of generic time-independent nonlinearities [23, 0|) the 
STS will propagate without interactions. The purpouse 
of this paper is to provide a first systematic exploration 
of collisions of STS. A few results were already reported 
in Ref. . One of the main contributions of that paper 
was to realize that for a given set of parameters ajk it 
is possible to use STS to build explicit solutions of Eqs. 
0. These solutions are constructed by taking 

Uj ^ ajR^{r),j = 1, . . . ,n (6) 

for any set of coefficients aj satisfying 

ajial + ... + ajnof^ = 1, j = 1, . . . , n. (7) 

These solutions of Eqs. |(SJ) are called Stabilized Vector 
Solitons (SVS) because they correspond to solutions with 
appreciable overlapping of the different components Uj. 
In Ref. the stability of these structures as well as the 
way they arise from collisions between stabilized solitons 
was studied. 

In this paper we present many more examples of colli- 
sions between STS and the associated phenomenology. 



III. EFFECTIVE-PARTICLE MODEL FOR 
COLLISIONS OF STABILIZED TOWNES 
SOLITONS 

A. Motivation 

Before describing the direct numerical simulations of 
Eqs. |(SJ) in detail and the many different phenomenolo- 
gies observed we first present an effective-particle model 
for collisions of STS. This model will give us a few hints 
on the expected dynamics of the system. The idea, as in 
many other problems in Physics, is to assume that during 
the collisions stabilized solitons behave as particles in the 
sense that can be described qualitatively by a bell-type 
ansatz with a few free parameters. 

This type of assumptions allows to reduce the dy- 
namics to a finite number degrees of freedom and re- 
ceives many different names depending on the field 
of application: method of collective coordinates, aver- 
aged Lagrangian description, time-dependent variational 
method, effective-particle method, etc. 
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In our case we take as trial Gaussian functions, which 
is one of the standard choices 



Aj exp 



{v~Vj) 



2< 



(8) 



The i-dependent parameters have the following mean- 
ing: Aj is the wavefunction amplitude; the co- 
ordinates of the centroid; Wjx,Wjy are the widths along 
the X and y axis; vjx, Vjy the initial velocities and f3jx, Pjy 
are phase factors which are required to obtain reliable 
results . Although Gaussians do not have the same 
asymptotic decay as STS, our choice simplifies the calcu- 
lations while leading to the same qualitative relations for 
the parameters and, as we will see below, the resulting 
equations provide an elegant and simple picture for the 
dynamics of the centroids of the STS. 



B. Equations for the effective- particle parameters 

Eqs. ((S)) can be derived by means of a variational for- 
malism from a Lagrangian density which can be written 
as a sum of the Lagrangians for the linear operators plus 
the nonlinear interaction in the following simple form: 



1 " 



(9a) 

(9b) 
(9c) 



The standard calculations of the method 23] consist of 
minimizing the trial functions from Eqs. ^ over the 
Lagrangian density given by Eqs. The final result is 
a set of second order ordinary differential equations for 
the evolution of the effective-particle parameters. In the 
case of the centroids we obtain: 



E 

n 

E 



di, 



k 



k=l ■> 
" dl,k 



(lOa) 
(10b) 



With j k these equations are in the form of Newton's 
second law with Ijk playing the role of a potential ruling 
the interaction between pairs of solitons according to the 
expression: 



ivk-vj) 



ajkNkg(t) e 



ky 



"kx 



(11) 



ky 



being = J \uk\^dxdy the square norm of the fc-th 
wavepacket. The corresponding widths along x and y 



are determined by the following equations 



k=l 



1 " 

jy = ZJT- ~ E 



jy 



fc=i 



dwj. 



(12a) 
(12b) 



Finally, some complementary relationships (first integrals 
of motion) for the velocities and phase coefficients are 
also obtained 



Xj 'IXjPjxj 



Wi 



2wi 



"3V 



2w 



(13a) 
(13b) 

(13c) 
(13d) 



From Eqs. (|13|) it is evident that the quantities Vjx,Vjy 
play the role of the initial velocities of the distribution 
centroids. 



C. Fast modulation approximation 

In this paper we take the modulation g[t) = go + 
gi cos ^It although the same qualitative results are ob- 
tained with other choices for q. To get STS the period 
of g{t) must be very short 0)E)E3 and the mean value 
< g > must satisfy the relation < g >< —Nc, where Nc 
is the critical square norm for blowup. The oscillations 
of the wavepacket widths, induced by g(t), are very fast 
compared with the dynamics of the centroids ruled by 
the potential from Eq. (|ll|l with an effective range of 
the order of the sizes of the wavepackets. Therefore, as a 
first approximation, the evolution of Xj and yj can be de- 
coupled from the oscillations of ojjx and ojjy. With these 
assumptions the interaction between different wavepack- 
ets will be determined by a constant nonlinearity with 

g =< g >= go- 
On the other hand, the parameters of the modula- 
tion can be chosen to minimize the amplitude of the 
wavepacket oscillations jl3|. An example can be seen in 
Fig. ^ where the results of direct numerical integrations 
of Eqs. lO in the one-component case n — 1 are shown. 
The peak amplitude and width of the solution present 
a small variation when suitable parameters are taken. 
Thus, as a first approximation we can consider wavepack- 
ets of constant circular section Wjx{t) — Wjy(t) — Wj. 
With all these considerations the potential between pairs 
of solitons is given by 



^0 ajkNkgo -r^.jw'^ 

Ijk — 9 6 J*-' 



(14) 



where vo^ 
centroids. 



wj + wf. and rjk is the distance between 
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FIG. 1: Evolution of the width W = (/(s^ + y'')\u\'y/^ and 
of the peak amplitude A — max |it| of the numerical solution 
of Eqs. for n = 1 with a nonlinear coefficient (a) g(t) — 
-27r+87rcos(40t), (b) g{t) = -27r + 9.57r cos(40t). Insets show 
details of the fast oscillations for the amplitude (left panel) 
and for the width (right panel). 



unstable closed orbit (ro = V^w). Finally ii vq > vj 
closed orbits are not possible because the wavepackets 
move too fast to be captured. 

Taking this simple picture in mind, we will perform in 
the next section a numerical exploration of the dynamics 
of stabilized solitons. Our main interest will be to check 
the existence of orbits and other dynamical processes. 

IV. NUMERICAL SIMULATIONS OF 
STABILIZED SOLITONS COLLISIONS 

In this section we study the collisions of STS in the 
framework of Eqs. ISJ. We take initial data of the form 

u,{r,0) = RoA\\r-r,\\)e^''^-, (17) 

and study their evolution by means of a pseudospectral 
Fourier scheme with time evolution of split-step type 1^. 
The scheme incorporates absorbing boundary conditions 
to get rid of the small amounts of radiation which are 
generated both by the STS and by the coUisional pro- 
cesses. 

Now, we proceed to describe the different behaviors 
observed. 



Therefore, the interaction between two stabilized soli- 
tons is governed, in a first approximation, by an attrac- 
tive force which only depends on the distance rjk between 
the centroids. That is, the motion of the centroids will 
be similar to planetary mechanics, with a "gravitational" 
potential Ijf. decaying as a Gaussian instead as 1/r. In 
fact, several common properties of central forces like con- 
servation of angular momentum and the reduction of the 
2-body problem to a single body motion in the center 
of mass frame are straightforwardly derived. Specifically 
the reduced particle moves under the action of an effec- 
tive potential given by 

y^^i^e-V-V^ (15) 

where go < and J is the conserved angular momentum. 

Using elementary techniques it is easy to obtain a 
necessary condition for the existence of circular orbits, 
namely: 

CA2 = (16) 

where C — 2ajkNk\ga\/TTj'^ is a constant of motion and 
A = r^/w^. If C < Eq. (^Hl has no solutions. 

If C = e^/4 there is only one (Ao = 2) and for C > 
e^/4 there are two solutions (A^ < Aq and A„ > Aq). 
This means that depending on the initial velocity of the 
wavepackets vq two closed orbits can exit. The most 
internal one (rg < ^/2w) is stable whereas the external 
one {tu > V^w) is unstable. If vq is increased up to a 
threshold value v / (escape velocity) there exists only one 



A. Fast collisions 

As it was shown in Ref. (1^] when "fast collisions" 
of STS take place the solitons emerge with only slight 
variations of the amplitudes and widths. These be- 
haviors can be accounted for by the effective-particle 
model proposed in Sec. Illll To see this we focus on 
the collisions of two stabilized solitions initially placed 
at ri — (xi,yi),r2 — (2:2,2/2) with initial velocities 
vi = ivix,viy),V2 = iv2x,V2y)- FoUowiug the results of 
Sec. mil and assuming trial Gaussian functions of equal 
size {wix = W2x = Wx, wiy — W2y ~ Wy) we obtain the 
equation for the effective-particle parameters 




where £x = X2 — xi and iy — 2/2 — j/i. Moreover we 
have the complementary relations (|13|l and the conser- 
vation law N{t) = 'n:\A\'^LLixOJy = 7r| A(0)pcJa;(0)wj^(0). 



FIG. 2: [Color online] (a) Surface plots of |wi|^ and \u2\^ for 
times from t = to t = 3.4 corresponding to the numeri- 
cal simulation of Eqs. Q with g{t) = —2ti + 87rcos(40t), 
n = (-6,-6),r2 = (6,-6),i;i = (5/2l/^ 5/21/^) and 
D2 = (—5/2^''^, 5/2^''^). (b) and (c) show the comparison of 
the evolution of the widths calculated numerically (solid line) 
and from the effective-particle model (dashed line) using Eqs. 
TO . 

The different terms in Eqs. (|18|l account for the phe- 
nomenology observed in "fast coUisions" . For example, 
they contain an asymmetric interaction (notice the differ- 
ences between Eqs. (|18e|) and (|18tll ) due to the fact that 
the width evolutions depend on the separation of both 
solitons along the corresponding axes. One example of 
this kind of interaction can be seen in Fig. |21 where we 
plot the results of numerical simulation of Eqs. lO with 
n = (-6,-6),r2 = (6,-6),-«i = (5/21/2,5/21/2)^^2 ^ 
(—5/21/^,5/21/^). In this figure we also compare the 
results from Eqs. jSJ with those obtained from Eqs. 
(|18|l . We can see that although there are quantitative 
differences the effective-particle model reproduces the ob- 
served dynamics and the qualitative behavior of the sys- 
tem. 



B. Collapsing orbits 

We have studied the evolution of two stabilized solitons 
which initially are placed at ri = (—2, 0) and r2 = (2, 0) 
with initial velocities Viy = —0.06 and V2y = 0.06 along 
the y-axis. In this situation we have observed that the 



FIG. 3: [Color online] Surface plots of \ui\^ (solid) and |u2|^ 
(transparent) for different times corresponding to the simula- 
tion of Eqs. © with g{t) = -27r + Stt cos(40t), ri = (-2,0), 
r2 = (2,0), viy = -0.06 and V2y = 0.06. 



solitons move inwards on a spiraling orbit. Therefore, 
the distance between the centroids of the two wavepack- 
ets decreases monotonically and the solitons join at the 
origin periodically. 

This behavior is shown in Fig. 13 where we plot the 
evolution of juij^ and |u2p for different times. It can 
be also observed that the solitons interact continuously 
and that a partial splittin g ta kes place, leading to the 
formation of a pair of SVS UJ. 

The effective-particle method cannot take into account 
the readjustment phenomenon between stabilized soli- 
tons observed in Fig. |31 Therefore, this simplified de- 
scription is not valid any more and can be used only as a 
first approximation to the problem. A more sophisticated 
model will be presented in Sec. 



C. Expanding orbits 

For initial velocities larger than in the previous case 
it is possible to overcome the collapsing character of the 
orbit. If the initial velocity is above a threshold value 
(which is analogous to a escape velocity) the stabilized 
solitons follow outward trajectories. In Fig. 0] and Fig. 
Elwe show the results obtained for two stabilized solitons 
initially placed at ri = (—2, 0) and r2 = (2, 0) with input 
velocities viy = —0.08 and V2y = 0.08 along the y-axis. 





FIG. 4: [Color online] Pseudocolor plots of for different 
times corresponding to the evolution of component ui from 
simulation of Eqs. @ with g{t) — —2n + 87rcos(40t), ri = 
(—2, 0) and viy — —0.08. The last picture corresponds to the 
superposition of several snapshots from t = to t = 200. 



Again we observe the formation of SVS. Since what it 
is plotted is the amplitude of only one component 
or |M2p ) the existence of two main spots in many sub- 
plots of Fig. 21 and Fig. [3 show, once more, that the in- 
teraction between stabilized solitons not only affects the 
trajectories of their centers, but also induces a readjust- 
ment of the distributions and the dynamical formation 
of SVS. We will try to describe this phenomenon in the 
next section. 



D. Wavepacket splitting witli deflection 

In this subsection we consider another interesting case: 
one stabilized soliton initially at rest and another one 
approaching to it. The situation is shown in Fig. Eland 
Fig. [Tlwhere we plot respectively the evolution of ui and 
U2- The first wavepacket starts at ri = (0,0) with zero 
initial velocity and the second one is initially placed at 
r2 = (3, —3) with initial velocity V2y — 0.3 along the 
y-axis. It can be appreciated how the interaction leads 
to the formation of SVS by means of the same splitting 
mechanism observed in previous figures. Pictures show 
that both wavepackets split simultaneously and part of 
ui is dragged by one half of U2 forming a vector soliton. 
The remaining vector soliton is deflected. 



FIG. 5: [Color online] Same as Fig. |3 for the evolution of 
component U2 with r2 — (2, 0) and V2y = 0.08. 

E. Three interacting solitons 

We have also studied the three-body problem which 
corresponds to three solitons in the corners of an equi- 
lateral triangle of side d. Therefore, by taking the ori- 
gin of coordinates in the barycenter the solitons are 
placed at ri ^ (0,rf/V3), ra = (-d/2, -d/2V^), 
rs — (d/2, — (i/2V3). The initial velocities are vix = 
— wsin(7r/2), viy — 0, V2x — — 'ysin(7r/2 -I- 27r/3), V2y = 
w cos(7r/2 -I- 27r/3), v^^ = — w sin(7r/2 -I- 47r/3) and v^y = 
i;cos(7r/2 + 47r/3) being v the input velocity. The nu- 
merical results for d = 4 and u = 0.12 are shown in Fig. 
|S1 where we see the evolution of one of the wavepackets. 
The initial velocities are such that the attraction between 
solitons is not enough to keep them bounded and the 
wavepackets move outwards. As in the previous cases 
the distributions split during evolution, the initial dis- 
tributions readjust due to the interaction with the other 
two components and three two-part vector solitons are 
generated. 

V. THEORETICAL ANALYSIS OF THE 
WAVEPACKET SPLITTING 

A. Motivation and ansatz 

From the previous section it is clear that the effective- 
particle model used in Sec. IIIII cannot explain the wave- 
packet splitting observed in the most of our numerical 
simulations. Therefore, in this section we develop an- 
other approach to capture the main characteristics of 
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FIG. 6: [Color online] Pseudocolor plots of \ui\^ for different 
times corresponding to the evolution of component ui from 
simulation of Eqs. @ with g{t) — —2n + 87rcos(40t), ri = 
(0, 0) and vi = 0. The last pseudocolor plot corresponds to 
the superposition of several snapshots from t = to t = 50. 
It is also shown as a surface plot in the bottom-right picture. 



such splitting. We will use again an effective Lagrangian 
technique but now we consider a two-mode model de- 
scribed by the following ansatzs 



Ul 



All exp 



+Ai2 exp 



y 



U2 = All exp 



^A22 exp 



2«^L 


2<J 


ix -f if 


y" 






{x-lf 




2t^i. 




{x -f if 




2^ii. 





+ 



(19a) 



+ 



(19b) 



This choice now allows readjustment of mass between two 
parts related to each component of the vector system, i.e. 
Ul is now composed of two parts which can be used to 
emulate the splitting mechanism and the formation of 
SVS. 

However, the full analysis of Eqs. H19|) by effective La- 
grangian methods is cumbersome and to achieve some 
results we will make several simplifications which can be 
justified by the observation of the simulations. In the 
first place, when the initial stabilized solitons ui ans U2 
are equal, the splitting mechanism is a symmetric pro- 
cess meaning that the growth rate of one part of wi 
is the same that the growth rate of the corresponding 
part of U2. Therefore we can consider An = A22 = a, 



FIG. 7: [Color online] Same as Fig. |S| for the evolution of 
component U2 with r2 — (3, —3) and V2y = 0.3. 




t=6Q 



t=90 




FIG. 8: [Color online] Surface plots of |u3|^ for different times 
corresponding to the evolution of component U3 from simula- 
tion of Eqs. © with ri = {0,d/VS), ^2 = (-d/2, -d/2V3), 
ra — (d/2, — d/2\/3) and initial velocities vix = — t; sin(7r/2), 
viy — 0, V2x ~ — i'sin(7r/2 + 27r/3), 1121, = wcos(7r/2 + 27r/3), 
V3x = — i;sin(7r/2 -|- 47r/3), vsy — vcos{-k/2 + 47r/3) being 
d = 4,v = 0.12 and g{t) = -2tt + 87rcos(40t). 
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Al2 = A21 = P, Wix = W2x, Wiy = W2y, Wlx = W2x, 

Wiy — W2y Secondly, we will assume that the sizes of ev- 
ery part are very similar and remain close to their mean 
values except for the fast oscillations. So we consider that 

Wlx = W2x = Wix = W2x = Wx and Wiy = W2y = Wly = 

W2y — Wy being Wx,Wy constants. Finally, to get some 
insight on the problem we consider the case where the 
distance between centroids 2i is approximately constant 
which can happen for instance during a fast switching 
process. Thus, the only free parameters are a and 



B. Model equations 

The Lagrangian L = J Cdxdy obtained from the den- 
sity Lagrangian C given by Eqs. © when the ansatzs 
are as in Eqs. ((T^ is 

L = TTWxWy ^i{aa* — a* a + Pl3* — P* l3)+ 

+ i{a(3* — a* (3 + (3a* — (3*dt)e ™^ -f 
1 1 



+ (i«r + i/3r) 

+ (a/3* + /3a*) 



2wl 2wl 

w:' 



2wt 2wl 



l + e 



+ [6|ap|/3|2 + 2(a2(/3*)2 + (a*)2^2)]g 
+ 2|an/3|2 + 4(|a|2 + \P\''){a(3* + a*/3)e'^ M 



where we have taken an = 012 = 021 = 0,22 = 1- 

The standard calculations yield to a system of or- 
dinary differential equations for a and /3, 



1 e 
I e 1 



K2\ ( a 
2 U2 Ai j 1/3 



(21) 



where 



2N{t) 

-e -|- 



TTWxWy 



A2(t) 



+ (|a|2 + |/3|2)(H-e ^-2e ^) 
'wl-2e . 1 . \2Nit) 



2wt 2wl 



(l«l' + |/3n(2e 



-2) 



TTWxWy 



(22a) 



(22b) 



and the function N{t) is the square norm of the ansatzs 



System can be written as 



where 



^2{t) = 



vi V2\la 
V2 vx) \ii 



K2{t)e 



Ai(i) 



2(1 --e ^) 
Ai{t)e ^ - K2{t) 



2(1 -e" 



(23) 



(24a) 



(24b) 



Using system H23(l and its complex conjugate and tak- 
ing into account that ^1,1^2 G K can be immediately 
proved that 



d 

dt 
d 



^(|ap + |/3n = 0, 



-(a/3*+a*/3) = 0, 
and therefore three invariants exist 

i|2 



aP* + a* (3 
Nit) 



90, 
91, 
^0, 



(25a) 
(25b) 



(26a) 
(26b) 
(26c) 



being gi, 92, Nq constants. 

System (|23|) can be solved numerically to find the evo- 
lution of the amplitudes a and (3. Nevertheless we can 
solve it analytically by considering that i>i and 1^2 are con- 
stants since the only dependence on t is given by the non- 
linear coefficient g{t) and, as we discussed in Sec. IIII CI 
the dynamics is determined by an averaged nonlinearity 
.9 = 5o- By writing a ^ aj^ + iaj, P ^ (3ji + iffj with a^, 
aj, Pfj, Pi gM. the solutions of system (I23II are obtained 
by using basic techniques from the theory of ordinary 
differential equations and the result is 





( 


Pr 


1 

= -M 


ai 


2 


\PiJ 


V 




where 



M = 



cos(i'i 
sin(z^i 



f-Mi M2 -M3 M4> 

Ml -M2 -M3 M4 

M2 Ml M4 M3 

\-M2 -Ml M4 M3/ 



(27) 



and 



Ml 
M2 
M3 
Ma 



a/(0) ' 
ai?(0) 
a/(0) - 



PM, 
- Pr{0), 

pm, 

-Pr{0), 



(28) 



(29a) 
(29b) 
(29c) 
(29d) 



N{t) ^\\u^\\l= ^w^Wy[\a\'' + |/?|2 + {aP* +a*P)e 



being a_R(0), a/(0), /3fl(0), /3/(0) the initial conditions. 
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50 100 ^ 150 200 

FIG. 9: Comparison of the evolution of the peak amplitude 
of the wavepacket iii corresponding to the simulation of Fig. 
^and Fig. |3 obtained by direct numerical simulation of Eqs. 
with the periodic oscillations from Eqs. 1311 predicted by 
the wavepacket splitting model. 

C. Validation 

Taking as initial conditions ai^(O) = Q:/(0) — 
Pj^{0) = (3j{0) = from Eq. ^ we obtain 

ctR = ao cos I'lt cos 1^2^, (30a) 

f^R = — sini/i<:sini/2t, (30b) 

a/ — aosinr^itcosj/2i, (30c) 

/?/ = ao cos 1^1 1 sin 1/2 (30d) 

and the amplitude evolutions are straightforwardly de- 
rived 

|a| = aol cosi/2i|, (31a) 
|/?| = ao|sinz/2t|. (31b) 

Eqs. H31|l mean that the splitting mechanism is oscil- 
latory and periodic as can be seen in Figs. I3I7I In Fig. 
El we compare the evolution of the peak amplitude of the 
wavepacket ui corresponding to the simulation of Fig. ^ 
and Fig. [Sl(for which £ = 2), obtained by direct numer- 
ical simulation of Eqs. lO with Eqs. (|31|l . The value 
of 1^2 is calculated according to Eq. H24b|) with a slight 
correction in the value of £ which is taking as = 1.88. 
This correction is due to the fact that numerical simula- 
tions of Fig. 0]and Fig. Oare made with Townes solitons 
intial data whereas the analysis of the wavepacket split- 
ting is made under the assumption of Gaussian initial 
data. Therefore, since Townes solitons decay at infinity 
as exp(— r) and Gaussians do it as exp(— r^), for a fixed 
separation 2£, the overlapping between Townes solitons 
is greater than between Gaussians. For these reason to 
compare both situations it is necessary to take a smaller 



value of £ in the case of Gaussian data. We see that 
there is a very good agreement between both curves up 
to around t = 70. After that numerical simulations show 
that solitons repel each other and the theoretical analy- 
sis is not valid any more, because we supposed constant 
separation between solitons. Therefore, we can conclude 
that our two-mode model is valid to predict the read- 
justment of mass between the wavepackets provided that 
the distance between parts remains nearly constant. As 
in the one dimensional case, the final repulsion of the 
wavepackets can be explained by taking into account that 
the relative phases between coherent solitons tend to sep- 
arate them |2J|. In this case, the solution of Eqs. l|5n|l 
gives a phase difference which takes the values ±7r/2. 

Finally, we must notice that the formation of vector 
solitons is the dominant tendency of the system in the 
slow collision regime. Therefore, the effective-particle 
model must be combined with the analysis of the par- 
tial splitting in order to get a more accurate picture of 
the dynamics. In fact this kind of solitons exhibit a tele- 
portation behavior as they suddenly vanish to appear in 
another place. This effect is very fast and can have po- 
tential applications in optical information processing. 



VI. CONCLUSIONS 

In this paper we have presented a detailed study of 
interactions between wavepackets which are stabilized 
against collapse by means of a modulation of the non- 
linear term, i.e. stabilized solitons. 

We have studied their dynamics by direct numerical 
simulations of the model equations. In this study, we 
have found that the system presents collapsing and ex- 
panding orbits depending on the initial configuration as 
well as other phenomena as deflection of one wavepacket 
due to the attraction and split in several parts in the 
corresponding n-body interactions. 

We have also developed a theoretical explanation of 
the wavepacket splitting based on effective Lagrangian 
methods and corroborated that effective-particle models 
allow to obtain conclusions for this system only in very 
limited situations. 
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